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ABSTRACT 

Three structural finite element programs are compared with theory, experimental 
data, and each other to evaluate their usefulness for modeling the 
thermomechanical deflection of ion engine electrodes. Two programs, NAS TRAN and 
MARC, used a Cray XMP and the third, Algor, used an IBM compatible personal 
computer. The shape of the applied temperature gradient greatly affects off-axis 
displacement, implying that an accurate temperature distribution is required to 
analyze new designs. The use of bulk material constants to model the perforated 
electrodes was investigated. The stress and displacement predictions are shown to 
be sensitive to the temperature gradient and the Young’s modulus, and insensitive 
to number of nodes, above some minimum value, and the Poisson ratio used. The 
models are shown to be useful tools for evaluating designs. Experimental 
measurement of temperatures and displacements was identified as the most critical 
area for further work. 


NOMENCLATURE 

surface area 

1/2 the center to center hole spacing 
effective acceleration distance 
electron charge 
Young’s modulus 
adjusted Young’s modulus 
maximum thrust 

instantaneous ligament half-thickness 
mean ligament half-thickness 
total beam current 

global stiffness matrix 

distance along a ligament 

hexagonal side length from hole pattern 

atomic mass 

ligament shape factor 

generalized force vector 
outer radius of circular plate 
radius 

von Mises stress 
path length 
electrode thickness 
temperature 

nodal temperature vector 

nodal displacement vector 
net accelerating voltage 
axial displacement 
two dimensional 
three dimensional 


a coefficient of thermal expansion 

X Child's Law constant 
AT center to edge temperature change 
y thrust correction factor for multiply 

charged ions and beam divergence 
v Poisson ratio 

v* adjusted Poisson ratio 

0 a stress in the axial direction 

a c stress in the circumferential direction 

a r stress in the radial direction 

ACRONYMS 

CPU central processing unit 
EBC Digital Equipment Corporation 
IBM International Business Machines 

PC personal computer 

INTRODUCTION 


Ion engines can provide a specific impulse of 
3000 to 8000 sec and thrusts of 0.1 to 1 N for 
systems with tens of kW available per engine. 1 Ion 
engines are very fuel efficient for many orbital 
applications, especially orbital transfer maneuvers 
and planetary missions. 2 * 3 Orbit transfers and 
planetary missions require scaling ion engines from 
tens of kW to MW powers for Mars transfer 
propulsion. 

An ion engine accelerates charged particles 
via an electrostatic field that is maintained between 
two thin, perforated electrodes. The charged 
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particle exit velocity is determined by the 
accelerating voltage and the charge to mass ratio of 
the particle. The thrusters are normally operated 
near the current limit which is defined by Child's 
Law, 


xa s v; 


3/2 


Jb= 


d2 


( 1 ) 


where Jb is the total beam current and x, sometimes 
called Child’s Law constant, is primarily determined 
by geometry. 4 Converting the current to thrust 
gives, 


xyA ,v* 

Fmax= ~j2 



( 2 ) 


where y is a correction factor to account for doubly 
charged particles and beam divergence. 5 The 
detailed geometry of the electrodes, ie. hole size, 
hole alignment, electrode thickness, etc., modify the 
effective acceleration length, the Child’s constant 
and the thrust correction factor. The contributions 
have been studied and are described in detail in 
references 4 through 7. In addition to optimizing 
the detailed geometry, the thrust per engine may be 
increased without affecting the specific impulse by 
reducing the electrode gap to reduce the effective 
acceleration distance or by increasing the beam 
extraction area of the engine. 

If the gap is too small, arc discharges from the 
screen to the accelerator electrode may occur too 
frequently. When the electrode area is increased, 
the effect thermally induced motion has on the 
electrode gap must be examined. Methods 
considered for increasing the electrode area include 
scaling up current engine designs or redesigning 
with a non-circular axisymmetric shape 8 or a non- 
axisymmetric shape. 9 Initial attempts to scale from 
30 to 50 cm diameter electrodes resulted in beam 
current densities significantly less than expectations 
based on extrapolations from 30 cm diameter 
thruster performance. 10 Non-circular electrode 
shapes provide a means of reducing the span to gap 
ratio of the electrodes while increasing the area, but 
they also diverge from the designs which have large 
design and operational experience data bases. 


technologists. 11-13 Models of two and three 
electrode systems were developed. All three of the 
modeling efforts cited used experimental data to 
confirm portions of the models. None of the models 
addressed non-linear effects of geometry and 
material properties. 

The goal of modeling efforts at Lewis Research 
Center is to verify and expand current modeling 
methods. This includes: checking the codes against 

theoretical and experimental data, documenting and 
justifying the selection of mesh elements and 
geometry, and verifying the non-linear modeling 
capabilities of the software. Toward these ends, 
models were built on three different finite element 
codes to analyze the structure of the electrode 
assembly. Comparisons were made within each 
code to evaluate model structures. Models from 
different codes were compared with each other and 
with plate-shell theory analysis to evaluate 
accuracy. Experimental measurements were 
compared with code predictions to evaluate the 
limitations and effectiveness of the modeling. The 
ability of all three software packages to handle non- 
linear effects was verified through a literature 
search. 14-16 The non-linear capabilities are 
required for modeling of non-axisymmetric 
electrodes, plastic strains, temperature variable and 
non-isometric material properties, and 
manufacturing defects. The intent of this paper is 
to describe the method and results of using finite 
element programs to model thermally induced 
motion of ion engine electrode assemblies. 

APPARATUS 

Study of the structural response of ion engines 
involves models to predict stresses and 
deformations, and experimental evaluation of 
temperatures and deformations. Three 

commercially available software packages were 
used to model the ion engine. The major code 
characteristics are described herein and 
summarized in Table I. The temperature and 
displacement measurements were described in a 
previous report, 17 so only an overview will be given 
here. 

Code P m .ri pti .ttog 


Computer modeling provides the potential to 
quickly analyze new designs and build the data 
base necessary to implement them. 

Structural analyses of ion engines have been 
conducted by a number of ion propulsion 


Three different finite element, structural 

analysis codes were used to model thermally 
induced stress and deflection of ion engine 
electrodes. The codes were chosen because they 

were available and all support 
3D Cartesian, 2D axisymmetric, 
and non-linear material and 
geometry effects analyses. 

The general operation of 
the codes is similar. A finite 
element model of the structure 
to be analyzed is generated 
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using a pre-processor. Displacement boundary 
conditions are applied using a combination of the 
pre-processor and a text editor. Each code used a 
different method to apply the thermal data. The 
main processor then assembles a global stiffness 
matrix using the data, and solves Hooke's Law, 
modified to account for thermal strain, 

0=£(u-af). (3) 

Displacement, stress, and temperature predictions 
are generated and written to a file. Each code 
generates data differently. A post-processor is used 
to display and compare the results generated by the 
different codes. 

A general description of each codes 
capabilities and specific characteristics follows, the 
detailed modeling approach will be discussed in the 
next section. 

Algor This code is written for an IBM 
compatible PC with a math co-processor. 14 It is 
written in a modular fashion, with separate 
programs for pre- and post-processing, stiffness 
matrix generation, thermal analysis, buckling 
analysis, static stiffness matrix solution, and 
dynamic mass and stiffness solution. Each 
individual routine has a unique name. For 
simplicity, the collection of programs will be 
referred to as Algor. Temperature effects are 
applied through the pre-processor. Data is 
generated at the element nodes. The output files 
from the separate routines can be used as input for 
the related routines, eg. the output from a thermal 
analysis can be applied as part of the load condition 
for a structural analysis, allowing one to restart a 
simulation from almost any previous solution. Only 
the pre- and post-processor and the elements 
necessary for static structural analysis in 3D and 
axisymmetric 2D have been used to date. 
Additional post processing was done using several 
personal computer 

spreadsheet and graphing 
packages. 

MARC This code is 
written and supported by 
the MARC Analysis Research 
Corporation, and works on a 
variety of machines, 
including mainframes by 
Cray Research, DEC, and 
IBM. 18 The program 
supports thermal, buckling, 
static structural and dynamic 
structural analyses. 

Temperature conditions 
were applied at the solution 
time through a user written 
subroutine. Data are 

reported for integration 
points within each element. 

As with the Algor package, 


restart files can be written from any point in the 
analysis and data from thermal analyses can be 
used as input for stress analyses. Again, this 
feature was not tested. At LeRC, MARC is installed 
on a Cray XMP. Pre-processing was done entirely 
with PATRAN, 19 a separate program running on a 
VAX computer. Post-processing was done with 
PATRAN and several personal computer 
spreadsheet and graphing packages. 

N ASTRAN This code is widely known, in part 
because several versions exist. A version supported 
by the MacNeal-Schwendler Corporation was 
installed on a Cray XMP. 20 The program supports 
thermal, buckling, static structural, dynamic 
structural and limited fluid flow analyses. Only the 
static stress analysis solution set was used. Thermal 
conditions were applied to the nodes by editing the 
input file. Data are generated at the centroids of the 
elements. Pre- and post-processing were identical 
to that used for MARC, excepting the different 
report formats used by the two codes. 

Hardware Description 

Because some of the model predictions will be 
compared with ion engine data, a brief description 
of the ion extraction system as well as hardware 
used to measure temperatures and deflections is 
presented. The testing described was completed in 
1981 and is described more completely in reference 
17. 

The experimental tests were conducted using 
a J-series electrode assembly mounted on a 30 cm 
diameter divergent field ion thruster. The J-series 
assemblies were fabricated in the late 1970’s for 
planetary propulsion applications using mercury 
propellant at 3 kW power levels. 21 A typical 
thruster is shown in Fig. 1. The engine body is 26 
cm long and 34.6 cm in diameter, with a cathode 
protruding 4 cm into the chamber, along the 
centerline. The plasma 
discharge is the main source 
of thermal input to the 
electrodes. An estimated 10 

percent of the discharge 
power is radiated through 
the electrodes. 22 The two 
electrode, J-series assembly 
is shown in Fig. 2. It consists 
of six basic parts: a titanium 
mounting ring which 
attaches to the plasma 
discharge chamber, the 

positively charged screen 
electrode which faces the 
plasma, a negatively charged 
extraction or accelerator 

electrode, a set of insulators 

which provide electrical 
isolation between the screen 
and accelerator electrodes. 
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Molybdenum 


Screen Electrode 


Stainless Steel 
Collar 


Shadow Shields 


14.2 


14.4 


15.0 


16.6 


Titanium Mounting Ring 


Alumina 


18.3 

all dimensions in cm 


insulator 


(vented) 


open area fraction of 
the screen and 
accelerator 
electrodes were 0.67 
and 0.24, 

respectively. The 
insulator assembly 


consisted 

machine 

separated 

cylindrical 

alumina 


of two 
screws 
by a 
piece of 
which is 


FIGURE 2.— J-SERIES ELECTRODE ASSEMBLY CROSS SECTION. 


MOUNTING RING 


ACCELERATOR ELECTRODE 


SCREEN ELECTRODE 


FIGURE 3. — J-SERIES ELECTRODE ASSEMBLY COMPONENTS. 

and two electrode stiffening rings. The stiffening 
rings and electrodes were made of molybdenum. 

The mounting ring attaches to a bracket at twelve 
points. For laboratory tests, twelve screws on the 
bracket join the bracket to the discharge chamber 
and center the electrodes on the axis of the 
chamber. For the structural analysis, 
the force transmitted through the 
brackets was assumed to be negligible. 

This allows the structural model of the 
assembly to be decoupled from a 
thruster model. The screen electrode 
has a hexagonal pattern of 0.19 cm 
diameter ion extraction holes on 0.22 
cm centers (Fig. 3). The accelerator 
electrode is patterned similarly but 
with 0.11 cm diameter holes. The 


drilled to permit 
outgassing (Fig. 2). 
The shadow shields 
over the alumina 
prevent shorts due to 
backsputtered material, and carry 
no structural load. A collar around 
the screws fixes the height of the 
accelerator electrode’s stiffening 
ring and thus determines the 
initial electrode gap. Twelve 
insulators were positioned equally 
spaced around the edge of the 
accelerator electrode’s stiffening 
ring. 

Tests were conducted to 

measure the temperature 

distribution across the electrodes 
and mounting rings and the 
displacement of the center of the 

electrodes. A brief description of the apparatus 
follows and a more complete description can be 
found in reference 17. Thermocouples were 

mounted at 10 locations on the electrodes and 
mounting ring, as shown in Fig. 4. Temperatures 



FIGURE 5. - DISPLACEMENT MEASURING PROBE. 
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TABLE II: TEMPERATURE AND DISPLACEMENT OF A 
J-SERIES ELECTRODE ASSEMBLY. 

Th er mo- 

Thermo- 

Radial 

Temperatures 

, °C 

couple 

couple 

Position, 

Discharge Power, W 

Location 

Number 

c m 

210 

350 

450 


1 

0.0 

278 

353 

397 

Screen 

2 

5.0 

266 

346 

381 


3 

14.4 

194 

294 

324 


4 

15.8 

191 

291 

319 

Mounting 

5 

16.9 

99 

215 

239 

Ring 

6 

18.3 

83 

196 

216 


7 

15.8 

124 

222 

246 


8 

14.4 

137 

234 

259 

Accelerator 

9 

5.0 

189 

276 

310 


10 

0.0 

209 

292 

327 

Centerline Gap 

Reduction 

, cm 

0.0025 

0.0051 

0.0076 

Centerline Screen Motion 

i, cm 

0.030 

0.046 

0.053 

Centerline Accelerator Motion, cm 

0.028 

0.041 

0.046 

• Thermocouple numbers refer to Fig. 4 

• No ion beam extraction. 

• Initial electrode gap was 0.051 cm. 

• Data taken from Ref. 17. 
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FIGURE 6 . - SCREEN ELECTRODE MEASURING 
SYSTEM SCHEMATIC 


codes and predict the behavior of the 
electrodes. The simplest model was compared 
to analytical solutions obtained using plate-shell 
theory. 23 * 24 * 25 Elements were then added to the 
model incrementally to observe their effects, and 
develop confidence in the model and methods used. 
Finally, the full mounting ring assembly was 
Measurements are listed in Table’ II. "The” probes modeled and the predictions obtained by applying 
were mounted on a bracket downstream of the measured temperature gradients to the model were 
electrodes (Fig. 6). Contact between a probe and the 
accelerator electrode or a probe and a post on the 


were recorded as the electrodes heated. Two 
stepper motor driven probes were used to measure 
the displacement of the electrodes (Fig. 5). 
Equilibrium temperature and displacement 


screen electrode completed an electric circuit. The 
motion of the electrodes was determined by the 
motion of the probe necessary to make and break 
contact with the electrodes. 

MODELING APPROACH 

All of the structural analysis was conducted 
by applying temperatures, and using boundary 
conditions to confine the motion of the edge of the 
electrode to the original plane, and confined motion 
of the center node to the axis of symmetry. The 
boundary conditions did not generate any load, so 
all forces and displacements were due thermal 
expansion and the geometry of the assembly. For 
the 3D Algor models, a 90° section was modeled 
(Fig. 7), and boundary conditions were applied to 
enforce radial symmetry, ie. no bending was 
allowed around the radial edges, and the nodes in 
the cutting planes were confined to motion in those 

planes. 

Four 
model 

configurations 
were examined 
to test the 
finite element 


compared with measured displacements. 

Co nfigurations 



Particulars of each of the configurations are 
given below. In order to reference the models 

easily, the models have been given codes, the first 
part is a Roman numeral indicating the stage or 
configuration modeled, a letter indicates the 
program used — a for Algor, m for MARC, and n for 
NASTRAN, and the Arabic number indicates 
whether the model used a 2D axisymmetric or 3D 
Cartesian coordinate system (Table III). 

Flat Plate — Stage I A solid, flat, circular plate 
of molybdenum, 15 cm in radius, and 0.038 cm 
thick was the simplest model used (Fig. 8A). 
Analytical solutions are available from plate-shell 
theory for this geometry with radially varying 
temperatures loads. 18 ' 20 This configuration 
represented the most complex problem with 
published solutions that contained any elements of 
the thruster model. The problem was modeled 
using all three codes. It was modeled with Algor in 
2D axisymmetric coordinates and 3D Cartesian 
coordinates and was the only model used on 
NASTRAN. 

Flat Plate with stiffening ring — Stage II A 
molybdenum stiffening ring was added to the Stage 
I model, extending from 15 cm to 18 cm radius, 


FIGURE 7. - J-SERIES 90° 
SECTION MODEL. 


TABLE III: LIST OF MODELS AND SOFTWARE USED 


Flat 

Plate 

Flat Plate with 
Stiffening Ring 

Dished Plate with 
Stiffening Ring 

Full 

Assembly 

2D Algor 

Ia2 

1 1 a 2 

1 1 1 a 2 

not modeled 

3D Algor 

I a 3 

1 1 a 3 

1 1 1 a 3 

IVa3 

MARC 

Im3 

I Im 3 

I IIm3 

not modeled 

NASTRAN 

In3 

not modeled 

not modeled 

not modeled 
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A: STAGE I 
B: STAGE H 



FIGURE 8. - MODEL CROSS SECTIONS. 


0.318 cm thick (Fig. 8B). There were no published 
analytical solutions available for this geometry with 
a radially varying temperature load. The model 
was used with Algor in 2 and 3D and in 3D with 
MARC. 

Dished Plate with stiffening ring — Stage 
III.. With this model, the shape of the electrode is 
assumed to be a spherical surface section with a flat 
extension and a stiffening ring attached to the 
extension (Fig. 8C). The radius of curvature of the 
electrode is 48.2 cm which equates to a dish depth 
of 2.2 cm with a circle of 14.4 cm radius defining 
the edge. 21 The dimensions were altered slightly 
from the stage I and II models, a flat flange extends 
from 14.2 to 15 cm radius, and the stiffening ring 
from 15 to 16.6 cm. The dimensions describe the J- 
series screen electrode. The stiffening ring was 
0.318 cm thick and the remainder of the model was 
0.038 cm thick. The material properties were 
modified across the material inside a 14.2 cm radius 
to simulate different hole patterns, ie. no holes, 0.11 
cm holes and 0.19 cm holes. Equations derived by 
Horvay 26 were used to modify the Youngs modulus 
(E) and the Poisson ratio (v) to account for the hole 
patterns, 

Fh 

E*=(l-v*) c m (4) 


30 cm Assembly with dished electrodes — Stage 
IV.. The electrode assembly model was used only 
with Algor. The model builds on the Stage III 
model, using two dished electrodes connected by a 
circular, L-section mounting ring and insulators that 
were modeled with constant property plate 
elements (Fig. 8D). The dished electrodes used bulk 
material constants to simulate 0.19 cm diameter 
holes in the screen electrode and 0.11 cm diameter 
holes in the accelerator electrode. The mounting 
ring (Fig. 2) was modeled using material properties 
for titanium and a thickness of 0.10 cm. The plate 
elements used to model the insulators had a width 
of 10° at 18.3 cm radius and a thickness of 0.35 cm 
(Fig. 7). Using properties for cold rolled steel, the 
bending stiffness of the plate elements modeling the 
insulators was matched to the estimated stiffness of 
the spacing screws in the radial plane, necessarily 
making them stiffer in the circumferential plane. 
Except for the constraints on the circular edge, the 
same boundary conditions used with the other 
electrode models were applied to the stage IV 
model. Displacement constraints on the circular 
edge were applied only to the screen electrode. 
These boundary conditions are enough to eliminate 
rigid body motions without introducing additional 
stresses due to overconstraints, and results in all 
axial deformations being reported relative to the 
stiffening ring on the screen electrode. 

RESULTS 

The results of the computer models were 
checked for sensitivity to the meshing and input 
parameters, agreement between different programs, 
and continuity with the other models. The output 
from the stage I model was compared with 
analytical predictions from plate-shell theory. The 
predictions from the stage IV electrode assembly 
model were compared with measured 


v*=l- 


3+2n(l+v) 


h m — L 


L 

d/ 

h 


1 


( 5 ) 


( 6 ) 


L-0 


where E* is the effective Young’s modulus, v* is the 
effective Poisson ratio, 2h m is a mean ligament 
thickness, h is the instantaneous half-thickness of 
the ligament or web (Fig. 9), and n is a shape factor 
for the cross section of the ligaments, taken to be 1 
for a rectangular cross section. The model was used 
with Algor in 2 and 3D and in 3D with MARC. 

Unless otherwise stated, all of the data 
presented were computed with a Young’s modulus 
of 280 GPa and a Poisson ratio of 0.324, 
corresponding to solid molybdenum, across the 
entire electrode. 



FIGURE 9. - HOLE GEOMETRY AND 
DEFINITIONS. 
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displacements. The stage II and III models were 
transitional stages between the stage I and IV 
models. 

Input Sensitivity 

In finite element modeling, the level of detail 
incorporated in a model must be balanced against 
the numerical accuracy of the program, and the 
time required to generate results. The objective of 
this sensitivity analysis was to evaluate the 
meshing options and the number of nodes for 
accuracy and computational efficiency and to 
evaluate the effect of varying inputs like 
temperature gradients and material constants. 

Evaluation of meshing Each element in a 
model is defined by its nodes. The 2D axisymmetric 
elements each have two nodes, one at each end 
point. The 3D elements used have four nodes each, 
defining the comers of a quadrilateral. The number 
of nodes determines computer time required and, to 
a large extent, the accuracy of the model. The 
model must have enough nodes to approximate the 
shape of the stress and displacement distribution. 
However, as the node spacing decreases, truncation 
errors become more noticeable and the computer 
time required increases. 

The results of IIIa2 models with 5, 10, 20 and 
30 nodes across the electrode were compared. The 
number of elements in each 2D model was one less 
than the number of nodes. Displacement 
predictions using 10, 20 or 30 nodes are within 0.1% 
of each other, but 5 node predictions diverge along 
the centerline and at the edge (Fig. 10). Stress 

0.05 


<-> 0.04 

£ 

j§ 0.03 

o 

< 

2 0.02 

o 

J 

S 

< 0.01 
0.00 

0 5 10 15 20 

RADIUS, CM 

FIGURE 10. — VARIATION IN DISPLACEMENT 
PREDICTIONS WITH MODEL SIZE, IIa2 MODEL. 

predictions from the 10, 20 and 30 node models 
were also within 0.1% and the 5 node model 
predictions diverged along the centerline. 

With 3D models, appropriate meshing must be 
chosen in the radial direction and in the azimuthal 
direction. Runs with circumferential spacings of 
15°, 10°, 7.5°, 5° and 2.5° showed no difference in 
either stress or displacement predictions. A 5° 
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FIGURE 11. — DISPLACEMENT DUE TO IMPOSED 
TEMPERATURE GRADIENTS, STAGE III MODEL. 

spacing was chosen for all of the 3D models because 
it had the flexibility needed to model the insulators 
on the stage IV assembly model and required less 
than two hours to execute on an IBM PC. 

The 3D models were tested with 10, 15 and 20 
radial nodes on Algor and MARC. Close agreement 
was seen with the 2D axisymmetric predictions 
except along the centerline (Fig. 11). The 3D models 
consistently predicted less displacement of the 
center elements than the 2D models predict. The 
deviation is due to the formulation of triangular 
elements which are used in the center, and will be 
discussed in the plate-shell theory section. 

There is an additional choice to be made with 
the 3D models. The full assembly can be modeled 
with a 90° section (Fig. 7) by using boundary 
conditions to enforce the axial symmetry of the full 
model. The boundary conditions are simple when 
the X-Z and Y-Z planes are used to section the 
model, and good agreement with a full 360° model 
was found by Brophy 13 using Algor. Using the 90° 
section, the number of nodes can be reduced by 
nearly a factor of four. Because computer time 
increases with the square of the number of nodes, 
the time required can be reduced by about a factor 
of 15. The drawback is that only axisymmetric 
buckling can be modeled using a 90° section. The 
stage III, 90° section Algor models required about 
40 minutes to converge. The 360° MARC and 

NASTRAN models required less than one minute of 
CPU time on the Cray, so they were not modeled as 
90° sections. 

Temperature gradient sensitivity The only 
load applied to the model was a temperature 
distribution. The differential thermal expansion 
across each electrode caused stresses and 
deflections. Because the ultimate goal of the 
modeling was to predict the behavior of the 
electrodes without knowing the exact temperature 
distribution, one must know how the shape of the 
temperature function affects the prediction. 

Data were taken previously at ten locations on 
the electrode assembly. 17 Three locations on each 
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of the two electrodes and four locations on the 
mounting ring were sampled (Fig. 4). The Stage I to 
III models had a single expansion coefficient, so 
there could be no stresses caused by a bulk 
temperature increase. This situation allowed the 
temperature distribution to be simplified to a 
gradient across the electrode. The measured 
gradient across the electrode was 73°C at 
equilibrium with a 450 W discharge. However, 
during thruster start-up, temperature gradients 
over 170°C were measured across the inner 14.4 cm 
of the electrode. A 100°C gradient was chosen for 
use with all of the modeling. 

In order to thermally load the models, it was 
necessary to interpolate between the points to 
generate data at the coordinates the model 
required. Three different profiles were used to 
model the thermal load. The simplest temperature 
profile used was linear, ie 100°C at the electrode 
center to 0°C at 14.4 cm radius. The next profile 
was a quadratic curve that was fit to the three data 
points available for the screen electrode at 450 W 
equilibrium temperatures (Ref. 17). The fitted 
curve was then normalized to 100°C at the center 
and 0°C at 14.4 cm radius. Neither of these profiles 
reflect the physics of heat transfer, they merely fit 
the limited data available and the constraints 
chosen. The constants that were found are given in 
Table IV. 


TABLE IV: TEMPERATURE 

DISTRIBUTION EQUATIONS 

Equation Form 

A 

B 

C 

T(°C) = A + Br 

100 

-6.94 

i : 

TC°C) = A + Br + Cr 2 

100 

-3.02 

-0.272 


The third method of interpolating between the 
data available was used only with the stage IV, full 
assembly model. A cubic spline was fit to the data. 
The method and results are described below, in the 
section comparing the model results to the 
measured displacements. 

Both gradients from Table IV were applied to 
the stage three models. The results show less than 
3% variance for the displacement at the center, with 
measurable differences near the mid-radius. The 
quadratic distribution generated about 33% more 
displacement at the mid-radius than the linear 
profile caused (Fig. 11). 

Evaluation of material constant effects In 
order to model the perforations of the electrodes, 
bulk material constants were calculated for the 
modified molybdenum sheet. Equations 4 through 6 
were used 27 to calculate appropriate values of 
Young’s modulus and Poisson ratio. The equations 
define the new properties based on the hole 


TABLE V: 

MATERIAL PROPERTIES 

USED TO MODEL PERFORATED 

MOLYBDENUM SHEET. 

Geometry 

E* 

V* 

Q>en Area 


(Pa 


Fraction 

n o h oles 

280 

0.32 

0 

0.11 cm holes 

47 

0.33 

0.24 

0.19 cm holes 

8.1 

0.62 

0.67 



A. — DISPLACEMENT VARIATIONS. 



B. — STRESS VARIATIONS 
FIGURE 12. — SENSITIVITY TO POISSON RATIO, 
MODEL IIIa3. 

diameter and center to center spacing. Table V 
shows the new values for three geometries, no 
holes, 0.11 cm holes and 0.19 cm holes. As the hole 
diameter increases, Young’s modulus decreases and 
Poisson ratio increases. Physically, a Poisson ratio 
of greater than 0.5 is nonsensical for elastic material 
theory. Neither MARC nor NASTRAN will operate 
with the Poisson ratio defined for 0.19 cm holes. 
Using Algor and holding the Young’s modulus 
constant at 280 GPA, several cases were tried using 
the stage IIIa3 and IIIa2 models. Comparisons of 
the displacement predictions (Fig. 12a) show 
negligible variation for Poisson ratios representative 
of 0, 0.11, and 0.19 cm holes when the Young’s 
modulus is held constant at 280 GPa. Stress 
variations (Fig. 12b) are also negligible. The von 
Mises stress is a root mean square of the principal 
stress magnitudes that provides a simple means of 
comparing the stress vectors. 

There was also some question about the 
necessity of varying Young's modulus. Provided 
that there are no geometric or material 

nonlinearities in the structure, the Young's modulus 
should make little difference in the magnitude of 
the displacement. Indeed, for very simple 
structures, the Young’s modulus drops out of the 
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equation for displacement (see Eq. 10). But, when 
one adds materials with different thermal 
expansion rates, the system becomes nonlinear. 
Additionally, one of the goals of the modeling 

efforts is to simulate nonlinear effects due to plastic 
strain. This requires accurate stresses and more 
material information than is currently available. 
Holding the Poisson ratio constant at 0.324, the 

IIIa3 model was tested with Young’s moduli of 280, 
47 and 8.1 GPa for the perforated region of the 

electrodes. The prediction of the magnitude of the 
maximum stress is strongly affected by the value of 
Young's modulus chosen (Fig. 13). Because there 

was little difference in the radial motion of the 
nodes on the stiffening ring, there is little difference 
in the electrode displacement predictions. 
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FIGURE 13. — SENSITIVITY OF STRESS PREDICTIONS 
TO YOUNG’S MODULUS, MODEL IIIa3. 

When the Horvay modifications to the 
material constants were examined separately, stress 
and deflection were found to be undisturbed by 
changes in Poisson ratio, but Young’s modulus 
strongly affected the stresses. 

Comparison of Results with Plate-Shell 
Theory 

Thermal stress and deformation equations 
have been published for flat circular plates with 
arbitrary axisymmetric temperature gradients. 25 
The analyses describe stresses that occur when 
buckling is not allowed, and deformations for 
buckled cases. These analyses are compared to the 
predictions of the finite element programs for the 
same conditions. Comparisons of the stress 
predictions for an unbuckled flat plate allow a 
direct check on the stress predictions of the codes. 
The stress is directly coupled to the displacement 
through Hooke’s Law, and correct evaluation of the 
stress is important for future analyses involving 
non-linear material properties and plastic 
deformation. The displacement predictions were 
checked by comparing buckling mode shape 
equations 25 with the eigenvalues generated by the 
model. 


Without buckling For these analyses, a linear 
temperature gradient, with a maximum 

temperature of 100°C at the center and 0°C at 15 
cm, was used. This gradient causes stresses to 
develop due to the geometry of the plate. The 
equations allow only radial motion of the points 
within the plate. The material stresses generated 
under these conditions are: 


aEAT 

(7) 

0* ~ R)» 

aEAT 


o c = (2r - R), 

(8) 


and 


S= 3R V3r2 - 3rR + R2, (9) 


where a r and o c are the stresses in the radial and 
circumferential directions respectively, AT is the 
temperature difference between the center and 
edge, and S is the von Mises averaged stress. 23 For 
the following Stage I modeling, a Young’s modulus 
of 8.1GPa and a Poisson ratio of 0.324 were used. 

The stage I models are equivalent to the 
system modeled by Goodier 24 . The temperature 
difference was applied across the plate using all 
three finite element codes (Fig. 14). The Ia2 model 
agrees with the closed form theoretical analysis. All 
of the 3D analyses show errors at the center and 
smaller errors at the periphery. The error at the 
edge was deemed insignificant, but, because the 
region of most deflection and first likely contact is 
the center of the electrode, improvements to the 
center node stress predictions were sought. 



FIGURE 14. — COMPARISON OF FLAT PLATE STRESS 
PREDICTIONS WITHOUT BUCKLING, STAGE I MODEL. 

The reduced stress that the 3D codes predict 
at the center of the plate is a fault in the 
formulation of the triangular elements and is 
common to all 3D stress codes. 27 Several methods 
of reducing the effect of this error have been 
investigated. When the size of the triangular 
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FIGURE 15. - STRESS PREDICTION WITH SMALL 
CENTER ELEMENTS, MODEL lm3. 



FIGURE 16. — STRESS PREDICTION VARIATIONS WITH 
NODE DISTRIBUTION, MODEL Im3. 


elements was reduced by two orders of magnitude, 
the minimum stress and radius where the minimum 
stress occurs were drastically reduced (Fig. 15). 
These errors were likely due to truncation errors 
that are enhanced by the small dimensions, 

PA TRAN 19 provides a simple method of 
scaling the node distribution to produce elements 
that are graduated in size, with smaller elements 
near the center and larger ones at the edge. Models 
with four different rat,es of graduating the element 
sizes were tested using the Im3 model. 
Representative data from one of these is compared 
to the uniform mesh results, the theoretical 
predictions, and a model with 2.5 times the 
standard mesh density in Fig. 16. For all cases, the 
centerline stress predicted by MARC is lower than 
the theoretical predictions. Additionally, the 
minimum stress and point of minimum stress is 
incorrectly predicted when the element size is 
graduated. For the even, dense and graduated 

meshes, the stress predicted by the models was 
within 10% of the stress predicted with plate-shell 
theory. 

An early attempt at a solution was combining 
pairs of center triangles to make quadrilateral 


elements. The large angle that results along the 
circumference of the element leads to computational 
inaccuracies for the quadrilateral elements. These 
errors affect stress predictions past the mid radius, 
and destroy the azimuthal symmetry, making post 
processing more tedious. For these reasons, center 
quadrilaterals were deemed an unworkable 
solution. 

Over all, the accuracy of the models is good 
when neighboring elements are of comparable size. 
The model consists primarily of quadrilateral 
elements and, although the stress from the 

triangular elements affects the loading of the 
quadrilaterals, the stress predicted just two nodes 

away from the center is within 5% of the theoretical 
value. Any method of improving the accuracy of 
the models’ centerline stress predictions must 
eliminate dependence on the stress from the 

triangular elements. This can be done in the post 
processing by interpolating for the stress at the 

center and first radial nodes. Limited scaling of the 
elements and an increase in the number of elements 
reduce the influence of the center errors on the 
bulk predictions. Increasing the number of nodes 
also increases the number of points that can be used 
with the curve fit. 

With Buckling As with the analysis without 
buckling, a linear temperature gradient was applied. 
The equations model a flat plate that is simply 
supported at the edge. For these conditions, the 
mode shape is 


aAT 


(9r2R-4r3-5R3) ( 3^ } -(r^-R*) 0 ,'^ (10) 


where w is the displacement perpendicular to the 
plate 25 . The Poisson ratio must be less than 0.5 for 
this equation to hold. From the displacement 
function, the stress function can be calculated. The 
resulting stress is 


o Et 2 aATR 

S— O a — Q 


(R-r) 


( 11 ) 


where the only component is shear in the axial 
direction, a a , so the von Mises stress equals the 
shear stress. 

All three programs are capable of solving for 
the buckled mode shape. However only MARC was 
used to solve for the buckling mode. Comparison of 
the computed and analytical mode shapes (Fig. 17) 
shows reasonable agreement everywhere. The 
buckling solution is given by solving for the 
eigenvalues and eigenvectors of the generalized 
stiffness matrix. MARC normalizes the resulting 
eigenvector to one. Because the greatest deflection 
occurs at the center node, the center node solution 
must agree with theory. The fact that all of the 
other nodes also show good agreement with theory 
suggests that the solution at those points will not be 
strongly affected by the triangular element errors. 
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FIGURE 17. — COMPARISON OF PREDICTED 
BUCKLING MODE SHAPES, MODEL Im3. 


Transitional Configuration (II and III) 
Results 

Configurations II and III were used on Algor and 
MARC in order to build confidence in the similarity 
of the programs' output and link the analytically 
verifiable work to the experimentally verifiable in 
steps that permit an intuitive check on the model 
results. With the level II models, a stiffening ring 
was added to the flat plate. The stage III models 
incorporate the dished structure of the J-series 
electrodes and the stiffening ring from the stage II 
model. The models were loaded with a linear 
temperature gradient normalized to 100°C. Results 
are compared in Fig. 18. The stage II structure 
shows higher centerline stresses than the stage I 
structure and the minimum stress on the electrode 
occurs at a larger radius. The shape of the stress 
distribution across the electrode is quadratic for 
both stage I and II. The stage II stress distribution 
departs from the quadratic form at 15 cm, where 
the stiffening ring begins. The stage III structure 
deflects axially without the use of the buckling 
solution. The deflection allows a large reduction in 
the stress throughout the model. The sharp changes 
in the stress values at 15 cm radius is due to a 
change in thickness. The stage II and III models 
produced no inexplicable displacements or stresses. 

Comparison of Model IVa3 with Engine Data 
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FIGURE 18. — COMPARISON OF STAGE I, II AND III 
RESULTS MODELING ELECTRODES PERFORATED WITH 
0.19 CM HOLES. 


they were tried first. The program would not 
execute when pins were used in the 3D structure, so 
the model was rebuilt with 2D plate elements 
simulating the insulators. The program executed , 
but it is more difficult and less accurate to model an 
axisymmetric insulator with a single plate shell 
element rather than a pin element. 

Simple, linear and quadratic temperature 
profiles (Table IV), applied to the Stage I to III 
electrode models, were applied to the Stage IV 
assembly model to provide a basis for comparison. 
However, the Stage IV model includes materials 
with different expansion coefficients and therefore, 
the possibility of stresses generated by constant 
temperature loads arises. The simple temperature 
distributions from Table IV do not model the 
temperature distribution across the mounting ring, 
nor do they model the bulk temperature increase. 
To handle these additional conditions, a cubic spline 
was fit to the temperature data. 17 The coordinates 
of the test data were converted to a path length 
along a cross section in order to reflect the path of 
thermal conduction. The path defined runs from 
the center of the screen electrode, through the 
mounting ring to the center of the accelerator 
electrode. A standard spline routine from the 
IMSL/Math subroutine library, CSAKM, 28 was 
chosen to minimize oscillations in the resulting 
curve. The equations generated are shown in Table 
VI. 


Algor was used for the 
full electrode assembly model 
because it allowed direct 
access to the stiffness matrix 
parameters through the pre- 
processor, simplifying the 
temperature loading. A choice 
had to be made in 
representing the insulators 
(Fig. 2). It seemed that one 
dimensional pin elements 
would be the best option, so 



TABLE VI. 

— SPLINE FIT TO TEMPERATURE DATA. 


FOR r 

FORs 


POSITION 

BETWEEN 

AND 

BETWEEN 

AND 

TEMPERATURE, °C 


0.0 

5.0 

0.0 

5.0 

397- 1.831s-0254s2-0.004s3 

SCREEN 

5.0 

14.4 

5.0 

14.6 

381-4.668(s-5.0)-0283(s-5.0)i+0.016(s-5.0)3 


14.4 

15.8 

14.6 

16.0 

324-5.755(s- 14.6)+4.843(s- 14.6)2-2345(s- 14.6)3 

MOUNTING 

15.8 

16.9 

16.0 

18.2 

319-5.984(s- 16.0)-37.832(s- 16.02+10.920(s- 16.0)3 


16.9 

183 

182 

21.7 

239- 13.891(s- 182)+1583(s-18.2)2-t0.145(s- 18.2)3 

RING 

183 

15.8 

21.7 

29.0 

216+2^29(s-21.7)-0.084<s-21.7)2+0.041(s-21.7)3 


15.8 

14.4 

29.0 

30.4 

246+7. 882( s- 29.0) +4.077( s- 29.0) 2 1 96( s- 29.0) 3 

ACOURATCR 

14.4 

5.0 

30.4 

40.4 

259+6.384( s- 30.4) - 0.089( s- 30.4) 7- 0.002s- 30.4)3 


5.0 

0.0 

40.0 

45.0 

3 1 0+4 .02 1 ( s- 40.0) - 0.057( s- 40.0) 2- 0.0 1 3( s- 40.0) 3 


1 1 
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FIGURE 19. — ELECTRODE MOTION WITH SIMPLE 


TEMPERATURE GRADIENTS, MODEL IVa3. 
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FIGURE 20. — DISPLACEMENTS WITH MEASURED 
TEMPERATURE GRADIENT, MODEL IVa3. 



FIGURE 21. — DISPLACED SHAPE COMPARISON, 
MODEL IVa3. 

The results obtained using the simple linear 
and quadratic temperature gradients (Table IV) 
show identical motion of the screen and accelerator 
electrodes, indicating that the stiffening rings 
completely mask the effect of the expanding, 


titanium mounting ring (Fig. 19). The displacement 
of the accelerator also shows perfect, axial 
symmetry, indicating that the stiffening rings also 
mask the effect of the non-axisymmetric loads 
generated at the insulators. 

The displacements resulting from the spline 
fit are shown in Fig. 20 for a cross section at 30°. 
Again, the results were perfectly axially symmetric, 
so it made no difference whether the cross section 
was taken through an insulator or not. All of the 
electrode nodes displaced in the downstream 
direction, in qualitative agreement with measured 
results. 17 The gap between the electrodes was 
slightly reduced at the center. The center is shown 
to move less than the mid-radial points, which may 
be somewhat misleading. Because the electrodes 
are initially dished more than 2 cm, there is no 
reversal of the electrode curvature. In fact, a scaled 
plot showing the initial and final states (Fig. 21) 
shows no discernable displacement. 

Although the triangular elements still affect 
the displacements, the results from applying simple 
gradients (Table IV) show the maximum 
displacement occuring at the center (Fig. 19). 
Therefore, one must conclude that the reduced 
central displacement that the spline loaded model 
shows is due to the temperature gradient applied 
and not to the triangular element errors that plague 
the analyses. This conclusion suggests two 

possibilities, the electrodes may deflect in such a 
fashion, or the temperature load that causes such 
deflection is in error. More detailed experimental 
measurements of temperature and displacement 
could prove conclusively which is true. 

CONCLUDING REMARKS 

Three finite element codes were used to 
analyze the thermo-mechanical response of ion 
engine electrodes to various input parameters. The 
codes were compared with each other to examine 
the sensitivity to meshing choices and the 
differences between 2D axisymmetric and 3D 
models. The effects of applying different 
temperature profiles to the models were examined. 
The sensitivity of the stress and displacement 
predictions to Young’s modulus and Poisson ratio 
changes was examined to determine what material 
constant adjustments are necessary to model 
perforated molybdenum. The model predictions 
were compared with plate-shell theory predictions 
to evaluate the accuracy of the 2D and 3D results. 
In order to develop an intuitive understanding and 
continuity between simple plate-shell structures 
and a full thruster model, models with single 
electrodes and stiffening rings were used. Finally, 
results from a three dimensional model of the 
electrode assembly were compared to measured 
displacements. 

Several conclusions can be drawn from the 
work completed. Models with sufficient node 
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density to provide consistent electrode stress and 
displacements were demonstrated on both 
mainframe and personal computers. 

The shape of the temperature gradient 
applied affects the off axis displacements 
measurably. This means that reasonably accurate 
temperature distributions must be developed, 
through modeling or measuring, in order to analyze 
new designs. 

Perforated sheet can be modeled as solid by 
changing the Young’s modulus, Poisson ratio and 
yield stress. Yield stress sensitivities were not 
examined because no data beyond yield were 
available. Accurate values of Young’s modulus were 
essential for modeling perforated material; both 
stress and displacement predictions were affected. 
Poisson’s ratio was shown to have negligible effect 
on stress and displacement profiles. If the models 
are to be expanded to include plastic strain and 
other non-linear, non-isometric effects, materials 
testing to determine the properties is necessary. 

The two dimensional, axisymmetric models 
agree well with plate-shell theory. The 3D models 
all have errors where triangular elements are used. 
The errors are due to constant strain formulation of 


the triangular elements and cannot easily be 
eliminated. In spite of errors from the triangular 
elements in the 3D models, the displacement and 
stress predictions from the quadrilateral elements 
are very close to closed form solutions for flat 
plates. 

Qualitatively, displacement predictions made 
by the models agree with measured displacements, 
showing both electrodes moving downstream with 
the screen electrode moving further than the 
accelerator electrode. The existing experimental 
database is insufficient for strong quantitative 
comparisons to be made. Structural finite element 
models of ion engine electrodes are in a state of 
readiness. Further development and useful output 
await detailed experimental measurements. 
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